#include <iostream>
#include <fstream>
#include <vector>
#include <sstream>
#include <cstdlib>
#include <climits>
#include <cfloat>
#include <mpi.h>
#include "omd/omdtool.h"
#include "drift.h"
// For 2D problems...

using std::vector;
using std::ifstream;
using std::ofstream;
using std::istringstream;
using std::cout;
using std::endl;

using namespace omd;
using namespace edrift;
// equation parser
DriftSolver* myptr;
char err_msg[1024];

//---------------------FUNCTIONS---------------------

// table enumeration starting at 1
double read_table(double tnum, double x) {
  int tabidx=(int)tnum;
  double ret;
  if(tabidx<0||tabidx>=(int)myptr->table.size()) {
    string msg("equ: read to undefined table");
    msg.append(as_string(tnum));
    throw msg.c_str();
  }
  ret=myptr->table[tabidx]->read(x);
  return ret;
}

double dread_table(double tnum, double x) {
  int tabidx=(int)tnum;
  if(tabidx<0||tabidx>=(int)myptr->table.size()) {
    string msg("equ: read to undefined table");
    msg.append(as_string(tnum));
    throw msg.c_str();
  }
  return myptr->table[tabidx]->dread(x);
}

double neig_rcell_value(double i, double j, double k){
  int ii=myptr->current_i+(int)i;
  int jj=myptr->current_j+(int)j;
  int kk=myptr->current_k+(int)k;
  int idx=myptr->index(ii,jj,kk);
  if(idx<0) {throw "out of range";}
  return myptr->GetVector()[idx];
}

double neig_rsrc_value(double i, double j, double k){
  int ii=myptr->current_i+(int)i;
  int jj=myptr->current_j+(int)j;
  int kk=myptr->current_k+(int)k;
  int idx=myptr->index(ii,jj,kk);
  if(idx<0) {throw "out of range";}
  return myptr->GetSourceVector()[idx];
}

double neig_rcell_value_outrange(double i, double j, double k, double outr){
  int ii=myptr->current_i+(int)i;
  int jj=myptr->current_j+(int)j;
  int kk=myptr->current_k+(int)k;
  
  if(ii<0||ii>=myptr->nx) return outr;
  if(jj<0||jj>=myptr->ny) return outr;
  if(kk<0||kk>=myptr->nz) return outr;
  
  int idx=myptr->index(ii,jj,kk);
  return myptr->GetVector()[idx];
}

double neig_rsrc_value_outrange(double i, double j, double k, double outr){
  int ii=myptr->current_i+(int)i;
  int jj=myptr->current_j+(int)j;
  int kk=myptr->current_k+(int)k;
//  std::cout<<"HIT:"<<i<<","<<j<<","<<k<<"\n";
//  std::cout<<"HIT:"<<myptr->current_i<<","<<myptr->current_j<<","<<myptr->current_k<<"\n";
  if(ii<0||ii>=myptr->nx) return outr;
  if(jj<0||jj>=myptr->ny) return outr;
  if(kk<0||kk>=myptr->nz) return outr;

  int idx=myptr->index(ii,jj,kk);
  return myptr->GetSourceVector()[idx];
}

double ext_rcell_value(double p, double i, double j, double k){
  int ii=myptr->current_i+(int)i;
  int jj=myptr->current_j+(int)j;
  int kk=myptr->current_k+(int)k;
  // fixme: range checking
  int linkidx=(int)p;
  if(linkidx>=myptr->link.size()) throw "can not find link to external array";
  int idx=myptr->index(ii,jj,kk);
  if(idx<0) {throw "out of range";}
  return myptr->link.at(linkidx)[idx];
}

double ext_rcell_value_outrange(double p, double i, double j, double k, double outr){
  int ii=myptr->current_i+(int)i;
  int jj=myptr->current_j+(int)j;
  int kk=myptr->current_k+(int)k;
  
  if(ii<0||ii>=myptr->nx) return outr;
  if(jj<0||jj>=myptr->ny) return outr;
  if(kk<0||kk>=myptr->nz) return outr;
  
  int linkidx=(int)p;
  if(linkidx>=myptr->link.size()) throw "can not find link to external array";
  int idx=myptr->index(ii,jj,kk);
  return myptr->link.at(linkidx)[idx];
}

double extern_cell(double p) {
  int linkidx=(int)p;
  if(linkidx>=myptr->link.size()) throw "can not find link to external array";
  int idx=myptr->index(myptr->current_i,myptr->current_j,myptr->current_k);
  if(idx<0) {throw "BUG: WRONG CURRENT INDEX! WHY???";}
  return myptr->link.at(linkidx)[idx];
}

int cnti=0;
double test_fun(double a) {
  std::cout<<"called:"<<cnti++<<"\n";
  return a; //myptr->step;
}

//-----------------------------------------------------

void DriftSolver::equation_setup() {

  istringstream ss(replace_char(str_equ,'=',' '));
  string en,sequ;
  
  while(ss.good()) {
    ss>>en>>sequ;
    mu::Parser* equ=new mu::Parser;
    equpar.push_back(equ);
    equname.push_back(en);
    equexpr.push_back(sequ);
  }
  
  for(int i=0;i<(int)equpar.size();i++) {
    
    equpar[i]->DefineFun("Read", read_table); // read table
    equpar[i]->DefineFun("DRead", dread_table); // read table derivative
    equpar[i]->DefineFun("RCell", neig_rcell_value); // read cell at relative position
    equpar[i]->DefineFun("RSource", neig_rsrc_value); // source cell relative
    equpar[i]->DefineFun("Ext", extern_cell); // access to external linked data
    equpar[i]->DefineFun("ExtRCell", ext_rcell_value); // external linked data relative to cell coordinate
    equpar[i]->DefineFun("RCellOR", neig_rcell_value_outrange); // allow outofrange
    equpar[i]->DefineFun("RSourceOR", neig_rsrc_value_outrange); // allow outofrange
    equpar[i]->DefineFun("ExtRCellOR", ext_rcell_value_outrange);
    
    equpar[i]->DefineFun("Test",test_fun); // test
    equpar[i]->DefineConst("Pi",M_PI);
    equpar[i]->DefineConst("kb",8.6173324E-5); // eV/K
    equpar[i]->DefineVar("LX",&var_lx);
    equpar[i]->DefineVar("LY",&var_ly);
    equpar[i]->DefineVar("LZ",&var_lz);
    equpar[i]->DefineVar("DT",&timestep);
    equpar[i]->DefineVar("DL",&var_dl);
    equpar[i]->DefineVar("x",&var_x);
    equpar[i]->DefineVar("y",&var_y);
    equpar[i]->DefineVar("z",&var_z);
    equpar[i]->DefineVar("t",&var_t);
    equpar[i]->DefineVar("cell",&var_cell);

    equpar[i]->DefineVar("ix",&var_ix);
    equpar[i]->DefineVar("iy",&var_iy);
    equpar[i]->DefineVar("iz",&var_iz);

    equpar[i]->DefineVar("LOX",&var_lox);
    equpar[i]->DefineVar("LOY",&var_loy);
    equpar[i]->DefineVar("LOZ",&var_loz);
    equpar[i]->DefineVar("HIX",&var_hix);
    equpar[i]->DefineVar("HIY",&var_hiy);
    equpar[i]->DefineVar("HIZ",&var_hiz);
    
  }
  
  AssignEquVar();
  AssignTableVar();
  
  for(int i=0;i<(int)equexpr.size();i++) { // compiled equations
    equpar[i]->SetExpr(equexpr[i]);
  }
  
  var_lox=var_loy=var_loz=0.0;
  var_hix=nx*spacestep;
  var_hiy=ny*spacestep;
  var_hiz=nz*spacestep;
  
  equ_ready=true;

}

//----------------------------------------------------------------
// on return str will be cleared if out of time range or no range

DriftSolver::DriftSolver(int _nx, int _ny, int _nz, string _symbol, string _map) {
  nx=_nx;
  ny=_ny;
  nz=_nz; // if nz<=2 then 2D simulation
  if(nz<=2) nz=1;
  
  dim=nx*ny*nz;
  var_x=var_y=var_z=var_t=0.0;
  var_lx=var_ly=var_lz=0.0;
  var_lox=var_loy=var_loz=0.0;
  var_hix=var_hiy=var_hiz=0.0;
  var_dl=var_cell=var_step=0.0;
  timecounter=0.0;
  current_i=current_j=current_k=0;
  serial_process=false;
  strict_courant=true;
  memset(err_msg,0,1024);
  
  has_cell=false;
  has_src=false;
  has_ext=false;
  has_con=false;
  has_cap=false;
  
  algo=NULL;
  map=NULL;
  bound=NULL;
  T=A=S=Ext=NULL;
  dE=dS=dExt=sumE=sumS=NULL;
  cap=con=NULL;
  ncoupler=0;
  step=0;
  timestep=1.0;
  spacestep=1.0;
  name="drift-solver";
  neqvar=0;
  equpar.push_back(new mu::Parser);
  equname.push_back("main");
  equexpr.push_back("0");
  mainequ=equpar[0];
  equ_ready=false;
  link_updated=true;
  MPI_Comm_rank(MPI_COMM_WORLD,&me);
  MPI_Comm_size(MPI_COMM_WORLD,&nproc);
  initiate_mem();
  Remap(_symbol.c_str(),_map.c_str());
  
}

void DriftSolver::CheckStrings() {
  if(str_cell.size()) has_cell=true;
  if(str_ext.size())  has_ext =true;
  if(str_src.size())  has_src =true;
  if(str_con.size())  has_con =true;
  if(str_cap.size())  has_cap =true;
}

void DriftSolver::check_expression(string& str, int& t0, int& t1, eqntype& isequ) {

  if(str[0]=='$') { // execute in serial: every processor assign complete array
    serial_process=true;
    str.erase(0,1);
  } else serial_process=false;

  if(str[0]=='{') {
    isequ=equ;
//    std::cout<<str;
//    exit(0);
  }
  else if(str[0]=='[') isequ=compile;
  else isequ=number;
  
//  istringstream sequ(remove_char(remove_char(replace_char(str,':',' '),'{'),'}') );
  string sstr=replace_char(str,':',' ');
  sstr=remove_char(sstr,'{');
  sstr=remove_char(sstr,'}');
  sstr=remove_char(sstr,'[');
  sstr=remove_char(sstr,']');
  
  istringstream sequ(sstr);
    
  sequ>>str; // str holds only the equation or number

  if(sequ>>t0) {
    if(!(sequ>>t1)) t1=t0;
  } else {
    t0=-INT_MAX;
    t1=INT_MAX;
  }
  
}

double DriftSolver::evaluate(int ix, int iy, int iz) {
  double retv;

  if(!link_updated) UpdateLink();
  
  try {
    myptr=this;
    current_i=ix;
    current_j=iy;
    current_k=iz;
    
    var_ix=(double)ix;
    var_iy=(double)iy;
    var_iz=(double)iz;
    var_x=var_ix*spacestep;
    var_y=var_iy*spacestep;
    var_z=var_iz*spacestep;
    var_step=(double)step;
    
    int idx=index(ix,iy,iz);
    if(idx<0) throw "BUG: index out of range @evaluate(int,int,int) WHY??";
    
    var_cell=A[idx];
    retv=mainequ->Eval();
      
    // @mac: must be carefull using brackets
      
  } catch(mu::Parser::exception_type &e){
    string msg("equation error: ");
    msg.append(e.GetMsg());
    throw msg.c_str();
  }
  
  return retv;
}

double DriftSolver::evaluate(string& str) {
  double retv=0.0;
  if(!link_updated) UpdateLink();
  
  try {
    // non sense to use the following variable here
    myptr=this;
    var_x=var_ix=0.0;
    var_y=var_iy=0.0;
    var_z=var_iz=0.0;
    var_cell=0.0;
    var_step=(double)step;

    equpar[0]->SetExpr(str);
    retv=equpar[0]->Eval();
  } catch(mu::Parser::exception_type &e){
    string msg("equation error: ");
    msg.append(e.GetMsg());
    throw msg.c_str();
  }
  return retv;
}

// add to all parser the same variable list
void DriftSolver::AddEquVar(string vname, double value) {
  try {
    if(neqvar>=(FD_MAX_VAR-1)) throw "qllowed number of variable exceeded";
    memset(eqvar[neqvar].name,0,32);
    vname.copy(eqvar[neqvar].name,31);
    eqvar[neqvar].value=value;
    for(int i=0;i<(int)equpar.size();i++)
      equpar[i]->DefineVar(eqvar[neqvar].name,&(eqvar[neqvar].value));
    neqvar++;
  } catch(...) {
    string msg("Error adding variable: ");
    msg.append(vname);
    throw msg.c_str();
  }
}

void DriftSolver::AssignEquVar() { // expr -> list of var value pairs
  istringstream ss(replace_char(str_var,'=',' '));
  string vn,vv;
  while(ss.good()) {
    if(!(ss>>vn>>vv)) break;
    
    bool found=false;
    
    for(int i=0;i<neqvar;i++) {
      if(vn==eqvar[i].name) {
        eqvar[i].value=evaluate(vv);
        found=true;
        break;
      }
    }
    
    if(not found) AddEquVar(vn,evaluate(vv));
  
  }
}

void DriftSolver::initiate_mem() {
  
  // WARNING! dE is initiated on demand
  T=new double[dim];
  A=new double[dim];
  S=new double[dim];
  Ext=new double[dim];
  dS=new double[dim];
  dExt=new double[dim];
  con=new double[dim];
  cap=new double[dim];
  map = new int[dim];

  for(int i=0;i<dim;i++) {
    T[i]=A[i]=S[i]=Ext[i]=0.0;
    con[i]=cap[i]=1.0;
    map[i]=0; // all empty@init
  }
  
  counts=new int[nproc];
  displ=new int[nproc];
  recvnum=new int[nproc];
  sendnum=new int[nproc];
  lim=new int[2*nproc]; // lo&hi limits
  clist=new int*[nproc];
  nlist=new int*[nproc];
  glist=new int*[nproc];
  ndata=new double*[nproc];
  gdata=new double*[nproc];
  
  for(int p=0;p<nproc;p++) {
    clist[p]=NULL;
    nlist[p]=NULL;
    glist[p]=NULL;
    ndata[p]=NULL;
    gdata[p]=NULL;
  }  
  
}

DriftSolver::~DriftSolver() {
  if(map) delete[] map;
  if(T) delete[] T;
  if(A) delete[] A;
  if(S) delete[] S;
  if(Ext) delete[] Ext;
  if(dE) delete[] dE;
  if(dS) delete[] dS;
  if(dExt) delete[] dExt;
  if(sumE) delete[] sumE;
  if(sumS) delete[] sumS;
  if(cap) delete[] cap;
  if(con) delete[] con;
  if(bound) delete[] bound;
  
  for (int i=0;i<(int)table.size();i++) delete table[i];
  for (int i=0;i<(int)equpar.size();i++) delete equpar[i];

  delete[] counts;
  delete[] displ;
  delete[] sendnum;
  delete[] recvnum;
  delete[] lim;
  
  if(nproc>1) {
    for (int p=0;p<nproc;p++) {
      if(ndata[p]){delete[] ndata[p];ndata[p]=NULL;}
      if(gdata[p]){delete[] gdata[p];gdata[p]=NULL;}
      if(nlist[p]){delete[] nlist[p];nlist[p]=NULL;}
      if(glist[p]){delete[] glist[p];glist[p]=NULL;}
      if(clist[p]){delete[] clist[p];clist[p]=NULL;}
    }
    
    delete[] ndata;
    delete[] gdata;
    delete[] nlist;
    delete[] glist;
    delete[] clist;
  }
  
  if(wfile.is_open()) wfile.close();
}

bool DriftSolver::mycell(int ix, int iy, int iz) {
  int idx=index(ix,iy,iz);
  if(idx<0) return false;
  if(serial_process) return true;
  // if(!map[idx]) return false; // skip empty cells
  return (idx>=lim[me*2] && idx<=lim[me*2+1]);
}

void DriftSolver::select_ncell(int cond, int idx, char** nmap) {
  if(cond) {
    for(int p=0;p<nproc;p++) {
      if(p==me) continue;
      if(idx>=lim[p*2] && idx<=lim[p*2+1]) {
        nmap[p][idx]=1;
      }
    }
  }
}

void DriftSolver::send_receive_list() {
  MPI_Request req[2*nproc];
  MPI_Status  stat[2*nproc];
  int rq=0;
              
  for(int p=0;p<nproc;p++) { // send to p
    if(p==me) continue;
    if(recvnum[p]) {
      MPI_Send_init(glist[p], recvnum[p], MPI_INT, p, 1, MPI_COMM_WORLD, &req[rq++]);
    }
  }
  
  // fetch index list to send
  for(int p=0;p<nproc;p++) { // receive from p
    if(p==me) continue;
    if(sendnum[p]) {
      MPI_Recv_init(nlist[p], sendnum[p], MPI_INT, p, 1, MPI_COMM_WORLD, &req[rq++]);
    }
  }
  
  MPI_Startall(rq,req);
  MPI_Waitall(rq,req,stat);
  for(int i=0;i<rq;i++) MPI_Request_free(&req[i]);
}

// assign -1 on neigh empty
// empty neighbour subject to von neumann condition
// in this case the temperature is isolated in the target

void DriftSolver::get_neigcell_index(int idx, int& l, int& r, int& t, int& b, int& n, int& f) {
  int i,j,k;

  get_cellcoord(idx,i,j,k);
  CellBoundary *cc=&(bound[map[idx]]);
  l=i-1; r=i+1; t=j-1; b=j+1; n=k-1; f=k+1;
  
  if(l<0) {
    if(forced_periodic&1) l=nx-1;
    else if(cc->l==FDCELL_PERIODIC) l=nx-1;
  }
  
  if(r>=nx) {
    if(forced_periodic&1) r=0; 
    else if(cc->r==FDCELL_PERIODIC) r=0;
  }
  if(t<0) {
    if(forced_periodic&2) t=ny-1; 
    else if(cc->t==FDCELL_PERIODIC) t=ny-1; 
  }
  if(b>=ny) {
    if(forced_periodic&2) b=0; 
    else if(cc->b==FDCELL_PERIODIC) b=0; 
  }
  
  if(nz!=1) {
    if(n<0) {
      if(forced_periodic&4) n=nz-1; 
      else if(cc->n==FDCELL_PERIODIC) n=nz-1; 
    }
    
    if(f>=nz){
      if(forced_periodic&4) f=0; 
      else if(cc->f==FDCELL_PERIODIC) f=0; 
    } 
  } else n=f=0;
  
  l=cc->l?index(l,j,k):-1;
  r=cc->r?index(r,j,k):-1;
  t=cc->t?index(i,t,k):-1;
  b=cc->b?index(i,b,k):-1;

  if(nz!=1) {
    n=cc->n?index(i,j,n):-1;
    f=cc->f?index(i,j,f):-1;
  }
  
}

void DriftSolver::initiate_comm() {
   
  displ[0]=0;
  occup=0; // # cells pre proc
  for(int i=0;i<dim;i++) if(map[i])occup++;
  ncp=occup/nproc;

  if(me==nproc-1) { // last proc take responsible for tail
    ncp=occup-ncp*(nproc-1);
  }
    
  MPI_Allgather(&ncp,1,MPI_INT,counts,1,MPI_INT,MPI_COMM_WORLD);
  
  for(int i=1;i<nproc;i++) displ[i]=counts[i-1]+displ[i-1];
  for(int i=0;i<nproc;i++) clist[i]=new int[counts[i]];
  
  int hilo[2];
  
  {
    int nc=0;
    int cc=0;
    int i;
    int start=0;
    for(i=0;i<me;i++) start+=counts[i];
    
    for(i=0;i<dim;i++) {
      if(map[i]) {        
        if((cc++)<start) continue;
        if(!nc) hilo[0]=i;
        clist[me][nc++]=i;
      }
      if(nc==ncp) break;
    }
    hilo[1]=i;
  }
  
  if(nproc==1) {
    lim[0]=hilo[0];
    lim[1]=hilo[1];
  }
  
  if(nproc>1) {
    
    {
      int tmp[occup];
      MPI_Allgatherv(clist[me],ncp,MPI_INT,tmp,counts,displ,MPI_INT,MPI_COMM_WORLD);

      int ii=0;
      for(int p=0;p<nproc;p++) {
        if(p==me) {ii+=ncp;continue;}
        for(int i=0;i<counts[p];i++) clist[p][i]=tmp[ii++];
      }
    }
    
    MPI_Allgather(hilo,2,MPI_INT,lim,2,MPI_INT,MPI_COMM_WORLD);

    char** nmap;
    nmap=new char*[nproc];
    for(int p=0;p<nproc;p++) {
      nmap[p]=new char[dim];
      for(int i=0;i<dim;i++) nmap[p][i]=0;
    }
    
    // check neighbor owner
    // avoid multiple-include
    for(int i=0;i<ncp;i++) {
      int idx=clist[me][i];      
      int l,r,t,b,n,f;
      CellBoundary *cc=&(bound[map[idx]]);
      get_neigcell_index(idx,l,r,t,b,n,f);

      select_ncell(l>=0,l,nmap);
      select_ncell(r>=0,r,nmap);
      select_ncell(t>=0,t,nmap);
      select_ncell(b>=0,b,nmap);
      
      if(nz!=1) {
        select_ncell(n>=0,n,nmap);
        select_ncell(f>=0,f,nmap);
      }
      
    }
    
    for(int p=0;p<nproc;p++) {
      recvnum[p]=0;
      sendnum[p]=0;
      for(int i=0;i<dim;i++)recvnum[p]+=(int)nmap[p][i];
    }
    
    for(int p=0;p<nproc;p++) {
      if(recvnum[p]==0||p==me) continue;
      glist[p]=new int[recvnum[p]];
      gdata[p]=new double[recvnum[p]];
      int n=0;
      for(int i=0;i<dim;i++)
        if(nmap[p][i]) glist[p][n++]=i;
    }
        
    // free...
    for(int p=0;p<nproc;p++) delete nmap[p];
    delete[] nmap;
    
    // ready and sync
    
    int allnum[nproc*nproc];
    MPI_Allgather(recvnum,nproc, MPI_INT, allnum, nproc, MPI_INT, MPI_COMM_WORLD);
    
    for(int p=0;p<nproc;p++) {
      if(p==me) continue;
      int* ptr=allnum+(p*nproc);
      sendnum[p]=ptr[me];
    }
    
    // allocate send and index buffer
    for(int p=0;p<nproc;p++) {
      if(!sendnum[p]) continue;
      nlist[p]=new int[sendnum[p]];
      ndata[p]=new double[sendnum[p]];
    }
    
    send_receive_list();
  }
}                       

// remaping must be accompanied by setting/rereading cap&con
// if needed also src
// if the non-time-ranged cap/con expression, this sould be 
// done automatic (sync_data).

void DriftSolver::Remap(const char* newsym, const char* newmap) {
  
  if(nproc>1) {
    for (int p=0;p<nproc;p++) {
      if(ndata[p]){delete[] ndata[p];ndata[p]=NULL;}
      if(gdata[p]){delete[] gdata[p];gdata[p]=NULL;}
      if(nlist[p]){delete[] nlist[p];nlist[p]=NULL;}
      if(glist[p]){delete[] glist[p];glist[p]=NULL;}
      if(clist[p]){delete[] clist[p];clist[p]=NULL;}
    }
  }
  
  // empty newsym or newmap -> leave untouched...
  if(newsym) {
    if(bound) delete[] bound;
    SetSymbolMap(newsym);
  }
  
  // WARNING! initiate only if map is changed...
  if(newmap) {
    SetMap(newmap);
    initiate_comm();
  }
  
}

void DriftSolver::SetTimestep(double dt){
  timestep=dt;
}
  
void DriftSolver::SetSpacestep(double h){
  spacestep=h;
  var_dl=h;
  var_lx=(double)nx*h;
  var_ly=(double)ny*h;
  h2inv=1.0/(h*h);
}

int DriftSolver::search_symbol(char c) {
  for(int i=0;i<nbound;i++) {
    if(bound[i].sym==c) return i;
  }
  return 0; // @head -> empty cell '.'
}

char DriftSolver::map_symbol(int num) {
  if(num<0) return '.';
  if(num>=nbound) return '.';
  return bound[num].sym;
}
 
// example:
//  a ++++++

void DriftSolver::SetSymbolMap(string str) {
  istringstream ss(str.c_str());
  vector<CellBoundary> vb;
  CellBoundary cc;

  cc.sym='.'; // dot is used for empty cell: head of symbol list
  cc.l=cc.r=cc.t=cc.b=cc.n=cc.f=FDCELL_EMPTY;
  vb.push_back(cc);
  
  while(ss.good()) {
    string tok, sb;

    if(ss>>tok>>sb) {
      if(sb.size()<4)
        throw "invalid map symbol entry";

      for(int i=0;i<(int)tok.size();i++) {
        cc.sym=tok[i];
        cc.l=cc.r=cc.t=cc.b=cc.n=cc.f=FDCELL_EMPTY;

        if(sb[0]=='+') cc.l=FDCELL_OCCUPIED;
        if(sb[0]=='p') cc.l=FDCELL_PERIODIC;
        if(sb[1]=='+') cc.r=FDCELL_OCCUPIED;
        if(sb[1]=='p') cc.r=FDCELL_PERIODIC;
        if(sb[2]=='+') cc.t=FDCELL_OCCUPIED;
        if(sb[2]=='p') cc.t=FDCELL_PERIODIC;
        if(sb[3]=='+') cc.b=FDCELL_OCCUPIED;
        if(sb[3]=='p') cc.b=FDCELL_PERIODIC;
        if(nz>1 && sb.size()==6) { // 3D
          if(sb[4]=='+') cc.n=FDCELL_OCCUPIED;
          if(sb[4]=='p') cc.n=FDCELL_PERIODIC;
          if(sb[5]=='+') cc.f=FDCELL_OCCUPIED;
          if(sb[5]=='p') cc.f=FDCELL_PERIODIC;
        }
        
        vb.push_back(cc);
      }
    }
  }
  
  nbound=(int)vb.size();
  bound=new CellBoundary[nbound];
  for(int i=0;i<nbound;i++) {
    bound[i]=vb[i];
  }
  
}

void DriftSolver::SetMap(const char* str) {
  for(int i=0;i<dim;i++) {
    if(str[i]==0x0) throw "insufficient number of character in map string";
    map[i]=search_symbol(str[i]);
  }  
}

void DriftSolver::EnableDelta() {
  if(dE) delete[] dE;
  dE=new double[dim];
}

void DriftSolver::EnableTally() {
  if(sumE) delete[] sumE;
  if(sumS) delete[] sumS;
  sumE=new double[dim];
  sumS=new double[dim];
  memset(sumE,0,dim*sizeof(double));
  memset(sumS,0,dim*sizeof(double));
}

void DriftSolver::KeepSource() {
  if(dS) delete[] dS;
  if(dExt) delete[] dExt;
  dS=new double[dim];
  dExt=new double[dim];
  memset(dS,0,dim*sizeof(double));
  memset(dExt,0,dim*sizeof(double));  
}

// one call per table
// three tokens per table: table@file scale const
// ore with renaming: table@file scale const name

// FIXME! use new form of TableReader constructor, if possible
void DriftSolver::AddTable(string str) {
  istringstream ss(str);
  while(ss.good()) {
    string stab,tname;
    if(!(ss>>stab)) break;
    int at=stab.find('@');
    
    if(at!=stab.npos){
      tname=stab.substr(0,at);
      stab.erase(0,at+1);
    } else {
      int dot=stab.find('.');
      tname=stab.substr(0,dot);
    }
    
    double scl,c;
    if(!(ss>>scl>>c)) scl=c=1.0;
    
    // TableReader will ignore tname if it is not OMD format
    // and assume the table is of raw format
    
    TableReader *tab=new TableReader(stab,tname,scl,c);
    table.push_back(tab);
    
    tabname.push_back(tname);
  }
}

void DriftSolver::AssignTableVar() {
  for (int i=0;i<(int)tabname.size();i++) AddEquVar(tabname[i],i);
}

/**
 communicates data at neighboring borders
 */

void DriftSolver::communicate(double* ptr) {
  if(nproc>1) {
    
    MPI_Request req[2*nproc];
    MPI_Status  stat[2*nproc];
    int rq=0;
    
    for(int p=0;p<nproc;p++) {
      if(!sendnum[p]) continue;
      for(int i=0;i<sendnum[p];i++) ndata[p][i]=ptr[nlist[p][i]];
      MPI_Send_init(ndata[p],sendnum[p],MPI_DOUBLE,p,1,MPI_COMM_WORLD,&req[rq++]);
    }
    
    for(int p=0;p<nproc;p++) {
      if(recvnum[p])
        MPI_Recv_init(gdata[p],recvnum[p],MPI_DOUBLE,p,1,MPI_COMM_WORLD,&req[rq++]);
    }
    
    MPI_Startall(rq,req);
    MPI_Waitall(rq,req,stat);
    
    for(int p=0;p<nproc;p++) {
      if(recvnum[p]) for(int i=0;i<recvnum[p];i++) ptr[glist[p][i]]=gdata[p][i];
    }

    for(int i=0;i<rq;i++) MPI_Request_free(&req[i]);
 
  }    
}

/**
 collects data and stor in every processor's array
 */

void DriftSolver::collect(double* ptr) {
  if(nproc>1) {
    double *tmp=new double[occup];
    double *dat=new double[ncp];
    
    for(int i=0;i<ncp;i++) dat[i]=ptr[clist[me][i]];
    
    MPI_Allgatherv(dat,ncp,MPI_DOUBLE,tmp,counts,displ,MPI_DOUBLE,MPI_COMM_WORLD);
    
    int ii=0;
    for(int p=0;p<nproc;p++) {
      if(p==me){ii+=ncp;continue;}
      for(int i=0;i<counts[p];i++) ptr[clist[p][i]]=tmp[ii++];
    }
    
    delete[] tmp;
    delete[] dat;
  }  
}

void DriftSolver::collect() {
  collect(A);
  collect(dS);
  collect(dExt);
  if(dE) collect(dE);
  if(sumE) {
    collect(sumE);
    collect(sumS);
  }
}

void DriftSolver::AssignArray(double* ptr, int ix, int iy, int iz, double val) {
  int idx=index(ix,iy,iz);
  if(idx<0) return;
  ptr[idx]=val;
}

bool once=true;

// some token in str is removed if already out of range
void DriftSolver::AssignArray(double* ptr, string& str) {
  if(str.empty()) return;
  
  if(once) {
    for(int i=0;i<(int)equexpr.size();i++) {
      equpar[i]->SetExpr(equexpr[i]);
    }
    once=false;
  }
  
  string tok,val;
  double dval=0.0;
  int t0,t1;
  eqntype isequ;
  istringstream ss(str.c_str());
  str=""; // keep rests of expression: accumulated
  
  while(ss.good()) {
    if(ss>>tok>>val) { // fixme: fetch all rest line!!
      check_expression(val,t0,t1,isequ);
      
      if(isequ==compile){
        // find a compiled equation
        int i;
        for(i=1;i<(int)equname.size();i++) {
          if(val==equname[i]) {
            mainequ=equpar[i];
            break;
          }
        }
        if(i==(int)equname.size()) {
          string msg("compiled equation not found:");
          msg.append(val);
          throw msg.c_str();
        }
      } else mainequ=equpar[0];
      
      if(step>=t0 && step<=t1) { // only if in time range
        
        // if only one time: skip
        if(t0!=t1) {
          
          if(serial_process) str.append(tok+" $");
          else str.append(tok+" ");
          
          if(isequ==equ)
            str.append("{"+val+":"+as_string(t0)+":"+as_string(t1)+"} ");
          else if(isequ==compile) {
            str.append("["+val+":"+as_string(t0)+":"+as_string(t1)+"] ");
            //std::cout<<"hit:"<<str<<"\n";
          }
          else
            str.append(val+":"+as_string(t0)+":"+as_string(t1)+" ");
        }
        
        if(isequ==equ) mainequ->SetExpr(val);
        else if(isequ==number) {
          // std::cout<<"get here: "<<tok<<":"<<val<<"\n";
          dval=as_double(val);
        }
        
        if(tok[0]=='@') {
          int p=val.find('=');
          istringstream scor(replace_char(val.substr(0,p),',',' '));
          val.erase(0,p+1);
          
          int ix,iy,iz,jx,jy,jz;        
          if(!(scor>>ix>>iy>>iz)) throw "wrong coordinate";
          if(scor>>jx>>jy>>jz) { // rectangle block
            // allow index out of range
            
            if(ix<0)ix=0;
            if(iy<0)iy=0;
            if(iz<0)iz=0;
            if(jx>=nx)jx=nx-1;
            if(jy>=ny)jy=ny-1;
            if(jz>=nz)jz=nz-1;
            
            if(nz==1){iz=0;jz=0;}
            
            // fill a block in range
            for(int i=ix;i<=jx;i++)
              for(int j=iy;j<=jy;j++)
                for(int k=iz;k<=jz;k++) 
                  if(mycell(i,j,k))
                    AssignArray(ptr,i,j,k,(isequ==number?dval:evaluate(i,j,k)));
          } else {
            if(mycell(ix,iy,iz)) // fill only one cell
              AssignArray(ptr,ix,iy,iz,(isequ==number?dval:evaluate(ix,iy,iz)));
          }
          
        } // if @ (range)
        else if(tok[0]=='*') { // fill all

//          if(isequ) mainequ->SetExpr(val);
//          else dval=as_double(val);
          
          for(int i=0;i<nx;i++)
            for(int j=0;j<ny;j++)
              for(int k=0;k<nz;k++) {
                if(mycell(i,j,k)) {
                  double v=(isequ==number?dval:evaluate(i,j,k));
                  AssignArray(ptr,i,j,k,v);
                }
              }
          
        } // if * (all cells)
        else { // symbol according to map

          // if(isequ) mainequ->SetExpr(val);
          // else dval=as_double(val);
          
          for(int i=0;i<nx;i++)
            for(int j=0;j<ny;j++)
              for(int k=0;k<nz;k++) {
                char c=map_symbol(map[index(i,j,k)]);
                if((tok.find(c)!=tok.npos) && mycell(i,j,k)){
                  AssignArray(ptr,i,j,k,(isequ==number?dval:evaluate(i,j,k)));
                }
              }
        }
      }
    }
  }
  
  // restore main equation
  mainequ=equpar[0];
}

void DriftSolver::sync_data() {
  
  int err=0,allerr;
  err_msg[0]=0x0;
  
  // clear arrays before to make them available after step
  // to fill source from program, use coupler->head()
  
  memset(S,0,dim*sizeof(double));
  memset(Ext,0,dim*sizeof(double));
  
  try {
    if(has_src)  AssignSource();
    if(has_cell) AssignCell();
    if(has_ext)  AssignExt();
    if(has_con)  AssignConductivity();
    if(has_cap)  AssignCapacity();
  
    for(int nc=0;nc<ncoupler;nc++) coupler[nc]->assign();
  } catch(const char* errst) {
    strncpy(err_msg,errst,1023);
    err=1;
  }
  
  MPI_Allreduce(&err,&allerr,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
  
  if(allerr==0) {
    communicate(con);
    communicate(cap);
    communicate(A);
  } else  throw err_msg;
}

bool dosync=true;

void DriftSolver::Inspect() {
  if(!equ_ready) equation_setup();
  sync_data();
  dosync=false;
}

// int dcnt=0;

void DriftSolver::Step(int n) {
  
  if(!equ_ready) equation_setup();
  
  // these are accumulated per call values
  memset(dS,0,dim*sizeof(double));
  memset(dExt,0,dim*sizeof(double));
  
  for(int loop=0;loop<n;loop++) {

    // maybe sync/inspected before loop (checking preloop values)
    // i.e. by calling Inspect()
    // this is allowed only one step per loop at the begining
    
    if(dosync) sync_data();
    dosync=true; // always sync while looping
    
    for(int nc=0;nc<ncoupler;nc++) coupler[nc]->head();
    
    if(algo) {
      memcpy(T,A,sizeof(double)*dim);
      algo->exec(); // NULL means no diffusion
    }
    
    // Source term enters here..
    // Value in source array is valid only for one step.
    // next step must be recreated
    // factor 2 is from averaging of conductivity
    //
    // Ext array is used for coupling purpose.

    for(int i=0;i<ncp;i++) {
      
      int idx=clist[me][i];
      if(!map[idx]) continue; // skip unmapped cells
      if(cap[idx]<=0.0) {
        string msg("Heat capacity may not be zero at this point (drift.cpp@1228)!! ");
        msg.append(as_string(idx)+" "+as_string(map[idx])+" "+as_string(cap[idx]));
        throw msg.c_str();
      }
      
      if(algo) {
        if(dE) {
          dE[idx]=(2.0*A[idx]*h2inv+S[idx]+Ext[idx])*timestep; // gives dE/vol
          A[idx]=T[idx]+dE[idx]/cap[idx];
        } else {
          A[idx]=T[idx]+(2.0*A[idx]*h2inv+S[idx]+Ext[idx])*timestep/cap[idx];
        }
      } else {
        if(dE) {
          dE[idx]=(S[idx]+Ext[idx])*timestep; // gives dE/vol
          A[idx]+=dE[idx]/cap[idx];
        } else {
          A[idx]+=(S[idx]+Ext[idx])*timestep/cap[idx];
        }
      }
      
      // in eV/A^3
      dS[idx]+=S[idx]*timestep;
      dExt[idx]+=Ext[idx]*timestep;
      
    }
//      std::cout<<"ed: "<<A[7501]<<"<---\n";
    
    // reserverd capability to call a coupler class
    for(int nc=0;nc<ncoupler;nc++) coupler[nc]->tail();
    
    var_t+=timestep;
    step++;
  } /* for loop */
  
  if(sumE) { // tally: if enabled
    for(int i=0;i<ncp;i++) {
      int idx=clist[me][i];
      sumS[idx]+=dS[idx];
      sumE[idx]+=dExt[idx];
    }
  }
}

double DriftSolver::GetSum(double *ptr) {
  double sum=0.0,ssum;
  
  for(int i=0;i<ncp;i++) {
    sum+=ptr[clist[me][i]];
  }

  MPI_Allreduce(&sum,&ssum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
  return ssum;
}

double DriftSolver::GetAverage(double *ptr) {
  return GetSum(ptr)/occup;
}

void DriftSolver::Dump(std::ostream &ofl) {
  if(me==0) {    
    ofl<<"# Drift Solver map file\n"
       <<"# (c) 2012, Yudi Rosandi\n"
       <<"# rosandi@gmail.com\n\n";
    
    ofl<<"dim "<<nx<<" "<<ny<<" "<<nz<<"\n"
       <<"timestep "<<timestep<<"\n"
       <<"spacestep "<<spacestep<<"\n"
       <<"symbol\n";
    
    for(int i=0;i<nbound;i++) {
            
      
      ofl<<bound[i].sym<<" ";         
      char l=bound[i].l==FDCELL_PERIODIC?'p':(bound[i].l==FDCELL_OCCUPIED?'+':'-');
      char r=bound[i].r==FDCELL_PERIODIC?'p':(bound[i].r==FDCELL_OCCUPIED?'+':'-');
      char t=bound[i].t==FDCELL_PERIODIC?'p':(bound[i].t==FDCELL_OCCUPIED?'+':'-');
      char b=bound[i].b==FDCELL_PERIODIC?'p':(bound[i].b==FDCELL_OCCUPIED?'+':'-');
      if(nz!=1) {      
        char n=bound[i].n==FDCELL_PERIODIC?'p':(bound[i].n==FDCELL_OCCUPIED?'+':'-');
        char f=bound[i].f==FDCELL_PERIODIC?'p':(bound[i].f==FDCELL_OCCUPIED?'+':'-');
        ofl<<l<<r<<t<<b<<n<<f<<"\n";
      } else  ofl<<t<<l<<r<<b<<"\n";
      
      /*
      ofl<<" # ("<<bound[i].l<<")"<<"("<<bound[i].r<<")"
      <<"("<<bound[i].t<<")"<<"("<<bound[i].b<<")"
      <<"("<<bound[i].n<<")"<<"("<<bound[i].f<<")";
      */
      
    }
    
    ofl<<"\n";
    
    int m=0;
    
    ofl<<"map "<<std::endl;
    for(int k=0;k<nz;k++) {
      for(int j=0;j<ny;j++) {
        for(int i=0;i<nx;i++) 
          ofl<<map_symbol(map[m++]);
        ofl<<std::endl;
      }
    }
    
    ofl<<"\n";
    ofl.flush();
  }

  MPI_Barrier(MPI_COMM_WORLD);

}

void DriftSolver::Dump(const string fname) {
  ofstream fl(fname.c_str());
  Dump(fl);
  fl.close();
}

void DriftSolver::CopyMap(char* cmap) {
  for (int i=0;i<dim;i++) {
    cmap[i]=map_symbol(map[i]);
  }
}

/**
 this function assign external link to a solver. This is called by external program (domain, console)
**/

void DriftSolver::Link(string pn, double* p) {
  link_updated=false; // mark it!
  // check if Link already exist
  for(int i=0;i<(int)linknm.size();i++) 
    if(pn==linknm[i]) {
      string msg("Link name already exists: ");
      msg.append(pn);
      throw msg.c_str();
    }
    
  link.push_back(p);
  linknm.push_back(pn);
}

/** 
 this function makes the index of link as variable name in equation.
 The domain will assign the name using "extern" directive. The access to the
 external array is done by Extern(link,i,j,k) function, where link is the 
 variable name, and i,j,k is the cell position.
 
 The call to this function must take place later, after all link is assign.
 Checked  @Step()
  
 **/

void DriftSolver::UpdateLink() {
  
  
  istringstream lss(str_link);
  
  while(lss.good()) {
    string pn;
    if(lss>>pn) {
      int i;
      for(i=0;i<(int)linknm.size();i++) {
        if(pn==linknm[i]) {
          AddEquVar(pn,i);
          break;
        }
      }
      
      if(i==(int)linknm.size()) {
        string msg("can not find external link name: ");
        msg.append(pn);
        throw msg.c_str();
      }
    }
  }
  link_updated=true;
}

// converting a real coordinate to cell index
int DriftSolver::coordtocell(double x, double y, double z) {
  int i=floor(x/spacestep);
  int j=floor(y/spacestep);
  int k=floor(z/spacestep);
  return index(i,j,k);
}

// get map from density array
// caller must take care that the array lengths are correct
// the length of 'syms' is used as length of 'limits'
// make sure that the symbols are listed in the model file!!!!

void DriftSolver::DensityMap(
                             double* dens, // density array. length=dim
                             const char* syms, // list of symbols, NULL terminated
                             double* limits, // limits of all symbols, sorted descending
                             const char* notouch // leave this character untouch or NULL
                             ) {
  
  int ns=strlen(syms);
  int nc=strlen(notouch);
  char cmap[dim];
  
  for(int i=0;i<dim;i++) {
    char ch='.',och=map_symbol(map[i]);
    
    for(int s=0;s<nc;s++)
      if(och==notouch[s]) {ch=och;break;}
    
    if(ch=='.'){
      for(int s=0;s<ns;s++)
        if(dens[i]>limits[s]) {ch=syms[s]; break;}
    }
    
    cmap[i]=ch;
  }
  
  Remap(NULL,cmap);

}

// convert map to array of double

void DriftSolver::MapToArray(double* dmap) {
  for (int i=0;i<dim;i++) {
    dmap[i]=(double)map[i];
  }
}

void DriftSolver::update_lox() {
  
  double tmp;
  bool getout;
  
  for(int i=0;i<nx;i++) {
    for(int j=0;j<ny;j++) {
      for(int k=0;k<nz;k++) 
        if(map[index(i,j,k)]) { 
          tmp=(double)i*spacestep;
          break;
        }
      if(getout)break;
    }
    if(getout)break;
  }
  
  MPI_Allreduce(&tmp,&var_lox,1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD);

}

void DriftSolver::update_loy() {
  
  double tmp;
  bool getout;
  
  for(int j=0;j<ny;j++) {
    for(int i=0;i<nx;i++) {
      for(int k=0;k<nz;k++) 
        if(map[index(i,j,k)]) {
          tmp=(double)j*spacestep;
          break;
        }
      break;
    }
    break;
  }
  
  MPI_Allreduce(&tmp,&var_loy,1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD);

}

void DriftSolver::update_loz() {
  
  double tmp;
  bool getout;
  
  for(int k=0;k<nz;k++) {
    for(int j=0;j<ny;j++) {
      for(int i=0;i<nx;i++) 
        if(map[index(i,j,k)]) {
          tmp=(double)k*spacestep;
          break;
        }
      break;
    }
    break;
  }
  
  MPI_Allreduce(&tmp,&var_loz,1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD);

}

void DriftSolver::update_hix() {
  
  double tmp;
  bool getout;
  
  for(int i=nx-1;i>=0;i--) {
    for(int j=0;j<ny;j++) {
      for(int k=0;k<nz;k++) 
        if(map[index(i,j,k)]) {
          tmp=(double)i*spacestep;
          break;
        }
      break;
    }
    break;
  }

  MPI_Allreduce(&tmp,&var_hix,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);

}

void DriftSolver::update_hiy() {
  
  double tmp;
  bool getout;

  for(int j=ny-1;j>=0;j--) {
    for(int i=0;i<nx;i++) {
      for(int k=0;k<nz;k++) 
        if(map[index(i,j,k)]) {
          tmp=(double)j*spacestep;
          break;
        }
      break;
    }
    break;
  }
  
  MPI_Allreduce(&tmp,&var_hiy,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);

}

void DriftSolver::update_hiz() {

  double tmp;
  bool getout;
  
  for(int k=nz-1;k>=0;k--) {
    for(int j=0;j<ny;j++) {
      for(int i=0;i<nx;i++) 
        if(map[index(i,j,k)]) {
          var_hiz=(double)k*spacestep;
          break;
        }
      break;
    }
    break;
  }
  
  MPI_Allreduce(&tmp,&var_hiz,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);

}

void DriftSolver::LinkCoupler(DriftCoupler* c) {
  ncoupler++;
  if(ncoupler>5) throw "Number of coupler over limit";
  c->init(this);
  coupler[ncoupler-1]=c;
  CheckStrings();
}

void DriftSolver::warning(string msg) {
  if(wfile.is_open()) {
    wfile<<msg<<std::endl;
    wfile.flush();
  } else std::cerr<<msg<<std::endl;
}

void DriftSolver::warning_file(string fname) {
  if(wfile.is_open()) wfile.close();
  wfile.open(fname.c_str());
}

// here calculate an optimum timestep
// by checking all cell against the Courant criteria
// better call this after a Step() call to make sure cap&con are actual

double DriftSolver::CourantAdapt(double factor, bool change_dt) {
  double crit=DBL_MAX;
  double newdt;
  
  for(int i=0;i<ncp;i++) {
    int idx=clist[me][i];
    if(!map[idx]) continue;
    if(con[idx]==0.0) continue; // this should not happen!!
    
    double cr=0.5*cap[idx]/con[idx]/h2inv; // dt<=h^2/2.D
    // std::cout<<cap[idx]<<"  "<<con[idx]<<cr<<"\n";
    if(crit>cr)crit=cr;
  }
  
  // here gather minimum!!

  MPI_Allreduce(&crit,&newdt,1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD);

  // std::cout<<"crit="<<newdt<<"\n";
  
  newdt*=factor;
  if(change_dt)timestep=newdt;
  return newdt;
}

void DriftSolver::SetAlgorithm(DriftAlgorithm* al){algo=al;algo->init(this);}
void DriftSolver::ClearAlgorithm(){if(algo) delete algo; algo=NULL;}


void DriftAlgorithm::exec() {
  
  double *A=solver->A;
  double *T=solver->T;
  double *con=solver->con;
  double *cap=solver->cap;
  double h2inv=solver->h2inv;
  int *map=solver->map;
  int me=solver->me;
  
  for(int i=0;i<solver->ncp;i++) {
    
    int idx=solver->clist[me][i];
    if(!map[idx]) continue;
    
    int l,r,t,b,n,f;
    solver->get_neigcell_index(idx,l,r,t,b,n,f);
    
    A[idx]=0.0;
    
    double K;
    
    if(con[idx]==0.0) continue;
    
    double crit=0.5*cap[idx]/con[idx]/h2inv; // dt<=h^2/2.D
    
    if(solver->timestep>crit) {
      static char st[1024];
      sprintf(st,"Courant stability criterion violated: "
              "step=%d time_step=%E crit=%E index=%d cap=%E con=%E h^2=%E",
              solver->step, solver->timestep, crit, idx, cap[idx], con[idx], 1./h2inv);
      // FIXME! strict_courant may belong to this class
      if(solver->strict_courant) throw st;
      solver->warning(st);
    }
    
    if(l>=0) if(map[l]) {
      K=con[l]+con[idx];
      if(K!=0.0) {
        K=con[l]*con[idx]/K;
        A[idx] +=K*(T[l]-T[idx]);
      }
    }
    
    if(r>=0) if(map[r]) {
      K=con[r]+con[idx];
      if(K!=0.0) {
        K=con[r]*con[idx]/K;
        A[idx] +=K*(T[r]-T[idx]);
      }
    }
    
    if(t>=0) if(map[t]) {
      K=con[t]+con[idx];
      if(K!=0.0) {
        K=con[t]*con[idx]/K;
        A[idx] +=K*(T[t]-T[idx]);
      }
    }
    
    if(b>=0) if(map[b]) {
      K=con[b]+con[idx];
      if(K!=0.0) {
        K=con[b]*con[idx]/K;
        A[idx] +=K*(T[b]-T[idx]);
      }
    }
    
    if(solver->nz!=1) {
      if(n>=0) if(map[n]) {
        K=con[n]+con[idx];
        if(K!=0.0) {
          K=con[n]*con[idx]/K;
          A[idx] +=K*(T[n]-T[idx]);
        }
      }
      if(f>=0) if(map[f]) {
        K=con[f]+con[idx];
        if(K!=0.0) {
          K=con[f]*con[idx]/K;
          A[idx] +=K*(T[f]-T[idx]);
        }
      }
    }
    
  } /*for i:ncp*/
}

